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1. Introduction. 

More than twenty years' experience of the lattice gauge community has taught us that it is 
always a good thing to have a bare action which respects symmetries, because then no fine tun- 
ings are required to preserve the symmetries at distances much greater than the cutoff scale. For 
this reason, essentially all lattice simulations of gauge theories perform simulations with gauge 
invariant lattice actions, and there is never any discussion about trading a small violation of gauge 
invariance in the simulation for larger volume or an apparently more efficient simulation algorithm. 
People have, however, always been willing to sacrifice chiral symmetry in their choice of lattice 
discretization. This seems somehow asymmetric: Why not perform simulations with lattice actions 
which preserve exact SU (Nf) x SU (Nf) chiral symmetry? 

The advantages of this approach are obvious: One does not have to separate the physical ex- 
plicit chiral symmetry breaking from a nonzero quark mass from the unphysical chiral symmetry 
breaking induced by lattice artifacts. The flavor content of the theory being simulated is unam- 
biguous. The index theorem is theoretically clean. The topological charge can be measured to be 
exactly what the dynamical fermions see during the simulation, not something which is determined 
by some post-processing procedure. And because the action preserves symmetries, correlation 
functions obey Ward identities which considerably simplify their theoretical analysis. For exam- 
ple, one does not have to spend any time measuring (and trying to remove) lattice-artifact additive 
mass renormalization or operator mixing. 

The way to do this is well known: use a lattice action which encodes the Ginsparg-Wilson[jl]] 
relation, an overlap [^, ^] action. This article is a summary of our experiences with simulations of 
two flavors of dynamical overlap fermions, using a version of the algorithm of Fodor, et al[Q]. It is 
a condensation of our two recent papers, Refs. ^ plus a little newer material. 

2. Algorithmic Issues 

Simulations with dynamical overlap have (at least) two problems: 

• They run so slowly 

• Changing topology is hard 

Our method of attack for the first problem is to replace the usual link variable gauge connection 
by a fat link. It has been known since at least 1998 [Q] that fat links improve the chiral properties 
of non chiral fermion actions (and the flavor symmetry properties of staggered fermions). The 
bottleneck has always been to find a smearing method which can be used in a molecular dynamics 
update, where the evaluation of the force requires a fat link which is differentiable with respect 
to its component thin links. Formulations like the Asqtad link solve this problem by "following 
the paths," but this does not give as much improvement as one would like. Our solution was 
provided by Peardon and Morningstar[||] with the "stout link" (invented in a Dublin public house): 
a multilevel blocking which is fully differentiable. In our runs, it pushes the thin link plaquette 
TrlJp ~ 1.7 up to about 2.8 (with two levels of smearing with p = 0.15). The number of Dirac 
operator matrix times vector multiplies per trajectory is reduced by about an order of magnitude 
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compared to simulations with a thin link gauge connection. The physics reason for this speedup 
is that fattening reduces the number of small eigenmodes of the kernel operator, improving its 
conditioning number. 

With a stout link, the fermions decouple from the UV fluctuations of the gauge field, and the 
mean size of the fermion force is reduced to about an order of magnitude smaller than the gauge 
force. This makes a multiple time scale integration algorithm very attractive. We run with the 
Sexton-Weingarten[^] form of this updating, taking the integration time step for the gauge fields to 
be 1/12 of that for the fermions. 

The square of the Hermitian overlap operator projected on one chiral sector is given by 



Hl{m)=2{Rl--^)Pa 



m 



2 



l + a££(A,-)|A,-)(A, 



Pa + m^ (2.1) 



with Rq the radius of the Ginsparg-Wilson circle, = i ( 1 + a/j ) the projector on chirality a and 
h{—Ro) the Hermitian kernel operator. The sign function e{h{—Ro)) is here given in its spectral 
representation. 

Because of the sign function in its definition, the effective action of the overlap operator has a 
discontinuity. It occurs when one eigenvalue of the kernel h{—Ro) changes sign during the molec- 
ular dynamics evolution. These are the surfaces in the space of the gauge fields on which the 
topology as defined by the index theorem changes by one unit. Ref. [Q] gives a prescription of how 
to account for this discontinuity in the HMC algorithm. One essentially measures the height AS of 
the step in the action (the potential of our Hamiltonian equations of motion) and if the momentum 
perpendicular to the surface is large enough one reduces it as one would do in classical mechanics. 
We will call this a "refraction" in the following. If the perpendicular momentum is too small, we 
flip it, and thus reflect the trajectory. With N the vector normal to the surface momenta n are thus 
updated by 



_ I -A^ {N\7i) sign(A^|7r) {N\tiY - 2ASf if{N\Kf>2ASf 
\-2N{N\n) if {N\7i)^ <2ASf 

The discontinuity AS of the effective action is caused by one eigenvalue changing sign, thus 
making the replacement 

Hl{m) -^hI± {4rI - m^)Pa\h) (Ao \Po (2.3) 
with |Ao) the zero mode. The corresponding step in the effective action can be evaluated using the 



Sherman-Morrison formula [10] 



(24) 

Interestingly, for the overlap not only can one compute the step in the effective action, but one can 
also give a closed form expression for the change in the fermionic determinant due to the change 
in topology: 

det^^(m) ,7 9^ , 1 , „ , 

;j^=l±(«J-„^){l<,|P,^P<,|A„>. (2.5, 
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In the actual simulation one faces the problem that the trajectories reflect most of the time off 
the zero eigenvalue surface and one never changes topology. The reason is that the change in the 
determinant from the starting configuration to the 'current' configuration in the MD evolution is 
only approximated well as long as the Dirac operator is similar to the starting one. The fluctuations 
are small as long as one has not changed topology. However, this is definitely not the case for 
the operator in a different topological sector. Since exp(— (^+l///^0) averages to the change in 
the determinant from the starting configuration to the end, large fluctuations mean that most of the 
time exp(— 0+l///^0) k, and the effective action is thus large, whereas only a few times do we 
get small effective actions. These two observations combine to give large \S most of the time, and 
thus a large auto-correlation time in the topological charge. 

To reduce fluctuations, we used the method proposed in Refs. [|ri|, [l^, which consists of 
rewriting the fermion determinant as as 



with mi = m and < with suitably chosen larger masses. In this method, only determinant 
ratios are evaluated using pseudo-fermions for the light quark masses. The change in the spectra 
while changing topological sector of the ratio H{m) /H{m') is expected to be less dramatic than the 
change of the spectrum of H{m). Only the determinant of H{mN) is evaluated directly. However, 
for a large mass mn the spectrum of is confined to a smaller region between mjj and 4Rq and 
the change in the spectrum therefore less drastic than for a smaller mass. One or two extra pseudo- 
fermion fields (Np = 2, 3 in Eq. |2.6| ) help some, but do not solve the problem. 

To quantify our difficulty, we compare in Fig. ^ the discontinuity of the effective action with 
the physical step from the fermion determinant. We subtracted the relevant quantity from the 
normal component of the momentum so that positive values correspond to reflections whereas the 
topology changes for negative values. We observe that the physical discontinuity would allow for 
significantly more changes in topology than the step in the effective action does. 

The low correlation between the estimator and the physical step height Eq. ( |2.5| ) shows up in 
the large auto-correlation time of the topological charge. Even though part of it is physics — lighter 
quarks make it harder to get from v=Otov = ±l — the height of the step grows with l/m^ instead 
of the expected determinant ratio, log m. Since the normal component of the momentum is roughly 
independent of the quark mass, it becomes more and more difficult to change topology. Indeed, 
Fig. ^ shows that the mean time between topological changes varies inversely with the square of the 
quark mass. The large auto-correlation time for the topology is a phenomenon that is also known 
with other fermions, e.g. improved staggered quarks. To the extent that these formulations know 
about topology, the step in the fermion action for the overlap might be replaced for them by a steep 
region which approximates the step. The result is the same: if the approximation of the determinant 
is bad, the step is overestimated most of the time and one does not change topology. 

3. The Topological Susceptibility 

In our second paper we made rough calculations of the topological susceptibility and chiral 
condensate using eigenmodes of the Dirac operator. We made simulations on 8^ lattices, at a lattice 



det//2(m) =det//2(mA,J Yl det 



(2.6) 



i=l 
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spacing of about a ~ 0.16 fm, with three quark masses, artiq = 0.03, 0.05, and 0.1. These lattices 
are really too small for physics, but they illustrate the useful features of a calculation with overlap 
fermions. 

We determined the string tension from the heavy quark potential and the Sommer parameter, 
from Wilson loops of temporal extent t = 2 and 3. The two measurements are not consistent, but 
we performed tests on 8'* and 12^^ quenched lattices which showed that the t = 3 potentials were 
consistent with ones from further separations. So we used their fit values. The lattice spacing varies 
by about ten percent as we change the quark mass. 

One picture. Fig. ^, illustrates our measurement of the topological susceptibility X- We take 
our measurements of ro/a and the topological charge time history to compute x^q- We have com- 
puted the lattice-to-MS' matching factor in perturbation theory and use it to convert the quark masses 



to their pi =1 GeV MS values. Durr[|13|] has presented a phenomenological interpolating for- 
mula for the mass dependence of the topological susceptibility, in terms of the condensate £ and 
quenched topological susceptibility Xq^ 

i = 4 + l (3.1) 

Taking E from our RMT analysis in the next section ( ro^Z ~ 0.43) produces the curve shown in 
the figure. 

Most published measurements of the topological susceptibility present them as a function of 
the pseudoscalar mass. Since we don't have spectroscopy, we can't do that. We can, however, use 
the Diirr curve as a fiducial, since most published measurements of the topological susceptibility 



present it, too. Our data (as well as that of Ref. 014|]) lies below the Diirr curve. Most measurements 



with nonchiral actions lie above it. (See, for example the figures in Ref. [||] or Since our 

quenched results give a value typical of simulations on larger lattices, % ~ (190 MeV)'*), we don't 
think we are seeing a finite volume effect. 
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Figure 1: The stochastic estimate of the height of the step compared to the actual change in the logarithm 
of the determinant from a subset of our ensemble. We subtracted the normal component of the momentum 
squared (which is typically less than 10) such that negative values mean refraction and positive ones reflec- 
tion. For mass ~ 0.03 on the left we have a number of events in the upper left quadrant that would have 
tunneled with the exact change of the determinant and only a few that actually tunneled (in the two lower 
quadrants). For nig = 0.05 the picture is similar, even though there are more tunneling events. 
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Figure 2: Monte Carlo simulation time between topology changes versus quark mass. 
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Figure 3: Topological susceptibility versus quark mass, in units of ro. The curved line is the Diirr interpo- 



lating formula, Eq. 3.1. The three horizontal lines give the quenched value and its error. 



4. The Condensate from Eigenmode Distributions 



It was proposed more than a decade ago that the distribution of the low-lying eigenvalues of the 
QCD Dirac operator in a finite volume can be predicted by random matrix theory (RMT) [|l^, [T8|, 
[T9| ]. Since then this hypothesis has received impressive support from lattice calculations, mainly 



quenched simulations pO], 21, 22, 23, 24, 23], but also some dynamical ones using staggered quarks 



6, 27] 



Typically, the predictions are made in the so-called epsilon regime, for which 1/A ^ L <C 
y/niji with A a typical hadronic scale. However, it has been found that they describe the data 
in a wider range. Two recent large scale studies, e.g., using the overlap operator on quenched 
configurations [24, 25], needed lattices with a length larger than 1.2 — 1.5 fm for RMT predictions 
match the result of the simulation. Our dynamical lattices have a spatial extent of about 1.3 fm. As 
we will see, random matrix theory describes our low-lying Dirac spectra quite well. 
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Our analysis is based on the distribution of the ^-th eigenmode from RMT as presented in 
Ref. [28] and successfully compared to simulation results in Ref. [^]. The prediction is for the 



distribution of the dimensionless quantity = pXk^V in each topological sector, with Xk the k- 
th eigenvalue of the Dirac operator, £ the chiral condensate, and V the volume of the box. The 
quantity p is the one-loop finite volume correction, p = 1 + c/{fjtLY where c is a "shape factor." 
These distribution are universal and depend only on the number of flavors, the topological charge 
and the dimensionless quantity niqiy . 

By comparing the distribution of the eigenmodes with the RMT prediction one can thus mea- 
sure the chiral condensate £. The main advantage of this method is that it gives the zero quark 
mass, infinite volume condensate directly. The validity of the approach can be verified comparing 
the shape of the distribution for the various modes and topological sectors. The main uncertainty 
comes from a too small volume which causes deviations in the shape, particularly for the higher 
modes. 

In Fig. we show the distribution of the two lowest eigenmodes of the overlap operator (scaled 
by rV) measured on the v = and v = ±1 parts of the arriq = 0.03 and aniq = 0.05 ensembles. 



We fit the RMT prediction from Ref. [ [28| ] to these distributions. The prediction agrees overall well 
with the measured distribution given the low statistics. However, the distribution of the lowest 
mode in the v = sector seems to have a tail at larger AZV that does not match the prediction. 
This could be an effect of the small volume. We also show the prediction for the distribution of the 
third mode from our fitted values of in the third column of Fig. ^ The RMT curve and the data, 
again, agree quite well. However for the | V | = 1 sector, the curve seems to be on the right of the 
data. This is probably a sign of the breakdown of RMT for eigenvalues larger than the Thouless 



energy [29, 



aniq 


pZrl 


0.05 


0.40(2) 


0.03 


0.44(2) 


0.01 


0.38(2) 



A combined fit to v = 0, 1, n = 1 , 2 at each mq (four distributions fit si- 
multaneously) gives the results shown in Table |l|. In our small volumes, and 
using the physical value for fj^ (93 MeV), p ~ 1.4, which is uncomfortably 
large. Dividing it out boldly gives Z ^ (280 MeV)^. 

After Ref. ^ appeared, we performed some simulations at lower quark 
Table 1: Condensate mass, aniq = 0.01. We restricted the topological sector to V = by forbid- 
^ ding refractions. It is unknown whether a particular topological sector is 

simply connected or not (in the latter case, forbidding topology changes 
might make the simulation non-ergodic). We ran two separate molecular dynamics streams to look 
for any obvious discrepancies. We did not see any: In the two streams, the plaquette and string 
tension parameters were consistent within statistical uncertainties; nothing looks unusual. We also 
did some running at am = 0.005 (or a 5 MeV quark mass), though not with enough statistics to fit 
the condensate. In both cases the code ran stably and quietly. 

Fig. ID shows the distributions and fit to the condensate from this data set. 

5. Conclusions 

Simulations with dynamical overlap fermions are poised to begin producing physics results. 
We are presently simulating on larger volumes in order to make a more reliable calculation of the 
condensate. We also continue to groom our algorithms. We believe that there are many more tricks 
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Figure 4: Distribution of the lowest two eigenmodes of the Dirac operator for our ensemble for v = 0,±1. 
The Unes are fits of the random matrix theory prediction to the data for the two lowest modes. The lines for 
the third mode are predictions. 



to be found and encourage others to work on dynamical simulations with overlap fermions. The 
physics payoffs are potentially very high. 

This work was supported by the US Department of Energy. 
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